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In an extensive computational experiment, we test Polyakov’s conjecture that under certain cir¬ 
cumstances an isotropic Heisenberg model can develop algebraic spin correlations. We demonstrate 
the emergence of a multi-spin U(l) order parameter in a Heisenberg antiferromagnet on interpene¬ 
trating honeycomb and triangular lattices. The correlations of this relative phase angle are observed 
to decay algebraically at intermediate temperatures in an extended critical phase. Using finite-size 
scaling, we show that both phase transitions are of the Berezinskii-Kosterlitz-Thouless type and at 
lower temperatures, we find long-range Ze order. 

PACS numbers: 75.10.-b, 75.10.Jm 


In statistical mechanics it is assumed [1, 2] that 2D 
Heisenberg magnets cannot develop algebraic order at 
finite temperatures since interaction of the Goldstone 
modes causes the spin-wave stiffness to renormalize to 
zero. However, in his pioneering work on this subject [3], 
Polyakov speculated that a 2D Heisenberg magnet might 
develop algebraic order if the system were to develop 
a “vacuum degeneracy”; he further suggested that this 
possibility might be explored experimentally. Recently 
Orth, Chandra, Coleman and Schmalian (OCCS) have 
proposed that frustration can provide a mechanism to re¬ 
alize Polyakov’s conjecture; here fluctuations induce an 
emergent XY order parameter that decouples from the 
rotational degrees of freedom [4, 5]. However these argu¬ 
ments were based on a long-wavelength renormalization 
group analysis, leaving open the possibility that short- 
wavelength fluctuations could preempt the scenario via 
unanticipated transitions into different phases [6-8]. In 
this Letter, we report a computational experiment that 
detects the development of an emergent XY order pa¬ 
rameter in a 2D Heisenberg spin model with power-law 
correlations, confirming the OCCS mechanism and its re¬ 
alization of the Polyakov conjecture. 

The OCCS mechanism relies on the formation of a 
multi-spin U(l) order parameter describing the relative 
orientation of the magnetization between a honeycomb 
and a triangular lattice. The development of discrete 
multi-spin order is well known in systems with compet¬ 
ing interactions: an example is the fluctuation-induced 
Z 2 order in the J\ — J 2 Heisenberg model [9]. This mech¬ 
anism is thought to be responsible for the high tempera¬ 
ture nematic phase observed in the iron-pnictides [10-13]. 
In the OCCS mechanism, the emergent 17(1) order pa¬ 
rameter is subject to a ’Lq order-by-disorder potential at 
short distances. At intermediate temperatures this po¬ 
tential is irrelevant (in the renormalization group sense) 


and scales to zero at long distances, leading to emer¬ 
gent power-law correlations. Remarkably, the stiffness 
of the emergent U(l) order parameter remains finite in 
the infinite system, despite the short-range correlations 
of the underlying Heisenberg spins. In this XY manifold 
the binding of logarithmically interacting defect vortices 
leads to multi-step ordering via two consecutive transi¬ 
tions in the Berezinskii-Kosterlitz-Thouless (BKT) uni- 



FIG. 1. (color online). Finite temperature phase diagram of 
classical windmill Heisenberg antiferromagnet as a function of 
inter-sublattice coupling Jth/J, J = VJttJhh- Below a copla- 
nar crossover temperature T cp , emergent XY spins appear and 
undergo two BKT phase transitions: at T> from a disordered 
to a critical phase with algebraic order and then at T< into 
a Z§ symmetry broken phase with discrete long-range order. 
At zero temperature the system undergoes a first order tran¬ 
sition at J t h = J from a 120° /Neel ordered windmill phase to 
a collinear phase. 
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versality class [4, 5, 14]. 

The Hamiltonian studied by OCCS is the “Windmill 
Heisenberg antiferromagnet”, given by H = H u + Hab + 
H tA + H tB with 

N 

H ab = j ab J2J2 s r s "+^’ (!) 

{-5a!,} 


where S" denote classical Heisenberg spins at Bravais 
lattice site j and basis site a £ {t, A, B}. The windmill 
lattice can be described as interpenetrating and coupled 
triangular (t) and honeycomb (A, B ) lattices. The indices 
5 ab relate nearest-neighbors of sublattices a, b , counting 
each bond once. The antiferromagnetic exchange cou¬ 
plings are J tt , J th = Jt A = JtB and J h h = Jab , and we 
introduce J = V Jtt Jhh- 

We employ large-scale parallel tempering classical 
Monte-Carlo simulations to obtain the finite tempera¬ 
ture phase diagram shown in Fig. 1 . As the emergent 
order parameter is a multi-spin object, we had to design 
a specific non-local Monte-Carlo updating sequence con¬ 
sisting of three sub-routines: (i) a heat bath step [15] in 
which a randomly chosen spin is aligned within the local 
exchange field of its neighbors according to a Boltzmann 
weight; (ii) a standard parallel tempering move [16, 17] 
for which we run parallel simulations at 40 temperature 
points and switch neighboring configurations according 
to the Metropolis rule; finally step (iii) is specifically tai¬ 
lored to our system where the emergent spins, defined 
below, exhibit a minute order-by-disorder potential. 
We select a (global) rotation axis perpendicular to the av¬ 
erage plane of the triangular spins, which exhibit (local) 
120 ° order, and rotate all honeycomb spins around this 
axis by a randomly chosen angle and accept according 
to the Metropolis rule. This Monte Carlo algorithm was 
applied at least for 9 x 10 5 Monte-Carlo steps of which 
the first half is discarded to account for thermalization. 

The emergent phases we are interested in occur for 
Jth < J where the zero temperature ground state is char¬ 
acterized by coplanar 120 ° order of the triangular spins 
and Neel order of the honeycomb spins (see Fig. 1) [18]. 
This order has SO(3) x 0(3)/0(2) symmetry and is de¬ 
scribed by five Euler angles {6,4>,%jJ) x (a,/3). As shown 
in the inset of Fig. 2, the angles (a, /3) describe the orien¬ 
tation of the honeycomb spins relative to the coordinate 
system i 7 (7 = 1,2,3) set by the triangular spins. The 
Euler angles ($,</>,?/>) relate i 7 to a fixed coordinate sys¬ 
tem. While the relative orientation can be changed with¬ 
out energy cost at T = 0, thermal fluctuations induce 
order-by-disorder potentials [19-21]. Considering Gaus¬ 
sian thermal fluctuations around the classical ground 
state, one finds a contribution to the free energy equal 
to [22, 23] 

= cos(2/3) [0.131 ^ - lO" 4 ^ cos 2 (3a) 


K 



T/J 


FIG. 2. (Color online) Coplanarity estimator k as a function 
of temperature for various values of Jth/J for system size 
L = 60. Inset shows definition of relative angles a and /3. 


The first term forces the spins to become coplanar 
(/3 = 7t/2) below a coplanarity crossover temperature 
T cp . More precisely, long-wavelength excitations out of 
the plane acquire a mass and are gapped out for T < T cp . 
The second term shows that the remaining U(l) relative 
angle a is subject to a Zg potential. 

As shown in Fig. 2, we track this coplanarity crossover 
within the Monte-Carlo simulations by measuring the 
coplanarity estimator 

3 J* 

k = i-— ^(cos 2 /3j), (3) 

j =1 


where cos / 3j = Sf ■ (S* x Sj +s± J with S tt being a nearest- 
neighbor vector on the triangular lattice. At high tem¬ 
peratures, where no relative spin configuration is pre¬ 
ferred, one expects k = 1/3, while for a completely copla¬ 
nar state holds k = 1. For uncorrelated 120° and Neel 
order at J t ,h = 0 one finds k = 0. Our Monte-Carlo re¬ 
sults show that coplanarity develops as soon as T < 0.25 J 
and k smoothly approaches unity for lower temperatures. 
Interestingly, n depends only weakly on J tb as long as 
Jth > J/ 10 . We define the location of the coplanar 
crossover T cp shown in Fig. 1 to be the location of the 
minimum of k. Note that down to the lowest tempera¬ 
tures we observe substantial out-of-the plane fluctuations 
and n < 1. We have identified these to be predominantly 
of short-wavelength nature. 

Below the coplanar crossover temperature T cp one may 
define emergent XY spins rrij at all Bravais lattice sites 
via projecting the honeycomb spin S j 4 (or S? ~ 
onto the plane that is spanned by the three nearest- 
neighbor triangular spins and normalizing 


(Sf-tuj-Sf-t 


c 2 ,j j 


(Sf ■ t ; j . Sf ■ Uj 


= ^cos OLj , sin Oij) 


( 2 ) 


(4) 
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We study the behavior of these emergent spins in the 
remainder of this paper. The local triangular triad t 7jJ is 
defined as follows: the spins on the triangular lattice are 
first partitioned into three classes {Sj’ X , Sj' Y , Sj' Z } as 

t X 

shown in Fig. 2. One then defines t-\ j = Sj and t 2 j to 

t y 

point along the component of <Sy that is perpendicular 
to tij. Finally, = tij xt 2 j completes the local triad. 
We show below that although the system exhibits out- 
of-the plane fluctuations and n < 1 , the emergent spins 
rrij decouple from these fluctuations and behave as U(l) 
degrees of freedom. 

To map out the low temperature phase diagram we 
analyze the correlations of the emergent spins rrij in the 
following. First, we define the total magnetization as 


m = 


1 

N 


N 

^2 rrij = |m|(cosa, sin a). 
j=i 


(5) 


The magnetization amplitude \m\ depends on the (lin¬ 
ear) system size L , in particular, it vanishes in the ab¬ 
sence of long-range order for L —► oo. Performing the 
Monte-Carlo average, we show the dependence of (|m|) 
with system size L in Fig. 3(a). While it vanishes faster 
than algebraic at large temperatures, it exhibits power- 
law scaling (|m|) oc L - ^ 7 ’)/ 2 with 0 < r/ < 0.3 for inter¬ 
mediate temperatures, a key signature of a critical phase. 
At the lowest temperatures, the exponent approaches 
zero and the magnetization saturates. To directly prove 
that the system develops (discrete) long-range order, we 
show the direction of the magnetization vector expressed 
as (cos( 6 a)) in Fig. 3(b). Clearly, (cos( 6 a)) approaches 
its saturation value of unity at low temperatures and 
large system sizes. The relative phase vector m points 
into one of the six directions preferred by the Xq potential 
in Eq. (2). The honeycomb spins are then aligned with 
one of the three triangular spin classes {S t,x , S tx , S'*’ 2 }, 
in agreement with the general order-from-disorder mech¬ 
anism that spins tend to align their fluctuation Weiss 
fields to maximize their coupling [ 21 ]. 

To determine the universality class of the phase tran¬ 
sition and the transition temperatures T> and T<, which 
partition the regimes of algebraic and long-range order¬ 
ing, we perform a finite-size scaling analysis of the XY 
susceptibility and magnetization for various values of 
Jth/J■ As shown in Fig. 4 we obtain perfect data col¬ 
lapse using a BKT scaling ansatz. Since the susceptibil¬ 
ity diverges when the system enters a critical phase, we 
can detect the upper transition at T> by analyzing 

x(T,T) = ^(|m| 2 ) = ^^m i f) (6) 

j 

for different temperatures T and system sizes L. We 
employ a BKT ansatz for the correlation length £> = 
exp(a>v / 7 ZA/i/T — T>) with a> being a non-universal 


(M) 
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FIG. 3. (color online), (a) XY magnetization amplitude (|m|) 
as a function of linear system size L for various temperatures 
T/J and fixed Jth/J = 0.8. On a double logarithmic plot 
it exhibits linear scaling within the critical phase with indi¬ 
cated floating exponent r)(T). It bends down in the disor¬ 
dered phase. Due to the finite system size we cannot clearly 
observe a saturation (at a finite value) at low temperatures, 
but ij approaches zero in a linear fit. (b) Direction of the 
magnetization expressed as (cos(6q)) as a function of T for 
Jth = 0.9J. A non-zero value signals breaking of the six-fold 
symmetry at low temperatures T < T<. Inset shows L = 12. 


constant [24-27]. Since \(T, oo) ~ £>(T) 2_r,> in the in¬ 
finite system, it holds that %(T, L) = L 2 ~ 71 Y X (^ > (T)/L) 
with a universal function Y x (x). For J t h = 0.6J we ex¬ 
tract the values T> = 0.200(4) J, a > = 1.9(3) and 7 y> = 
0.25(1) from optimizing the collapse. This agrees very 
well with the theoretically expected value 7 y> = 1/4 [14], 

Performing the analysis for other values of Jth yields 
data collapse of similar quality with a value r] > = 0.25 
within error bars. This determines T>(scal.) and the up¬ 
per phase transition line in Fig. 1 . As an independent 
way to determine T>, we use the power-law scaling of 
the magnetization with the system size L , which is ex¬ 
pected to be (|m|) oc L ~ v / 2 with 77 = 1/4 at the upper 
transition. This yields T>( 77 ) included in Fig. 1. The two 
temperatures agree within error bars with T > {rf) being 
systematically slightly larger. Finally, we note that we 
have also tried to achieve data collapse using a scaling 
ansatz corresponding to a second order phase transition, 
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but the resulting collapse is worse in this case, especially 
for data points close to the phase transition. 

To determine the lower transition temperature T< we 
perform a finite size scaling analysis of the magnetiza¬ 
tion amplitude (\m\(T,L)). Since it holds in the infi¬ 
nite system that (|m|(T)) oc £(T )<’ ,< ^ 2 with correlation 
length £< = exp (a < \/T < /y/T < — T ) and non-universal 
factor a<, it follows for a finite system that (|m|(T, L)) = 
L~’ n< / 2 Y rn {^ < (T)/L), where Y m (x) is a universal func¬ 
tion. In Fig. 4(b) we show the best data collapse for 
J t h = 0.6J which yields T< = 0.18(1), r] < = 0.11(1) 
and a< = 5.0(5). This is in good agreement with the 
theoretically expected value of r] < = 1/9 at the lower 
transition [ 6 , 14]. 

Two independent ways to obtain T< are (i) to in¬ 
vestigate the power-law scaling of (|m|) with system 
size and (ii) to directly look for the symmetry break¬ 
ing as indicated by the quantity (cos(6q;)). Using the 
first method, we find that our data can be fitted to 
log(|m|) a — ^^logL with a temperature-dependent 
slope r/(T) that is monotonically decreasing over the full 
range 0 < T < T>. At high temperatures, we find 
?y(T>) ~ 0.25 (as expected) and we define T<(? 7 ) as the 
temperature where 77 (T<) = 1/9. The fact that the sys¬ 
tem appears to be critical within our simulation even for 
lower temperatures (with an exponent r/ < 1/9) is a sim¬ 
ple consequence of the fact that the system size is much 
smaller than the correlation length [25]. If we were able 
to reach larger system sizes in the simulation, we would 
eventually see a saturation of (|m|) to a finite value. 

Next we discuss the second method to detect T<, 
namely direct observation of symmetry-breaking. We see 
in Fig. 3(b) that (cos( 6 a) approaches unity at low tem¬ 
peratures and large system sizes. In a finite-size system, 
we can observe this ordering only for not too small values 
of Jth > 0 . 8 J because the bare value of the order-from- 
disorder six-fold potential scales with ( J t h/J) 6 with an 
additional small numerical prefactor 10 -4 (see Eq. (2)). 
While the lower phase transition occurs when this po¬ 
tential becomes relevant at long lengthscales, indepen¬ 
dently of the bare value, the finite system size serves 
as a cut-off of the scaling making an effect of the po¬ 
tential only visible at sufficiently large bare values. To 
extract the transition temperature T< from (cos( 6 a)) we 
have to take into account that while at low temperatures 
the Gaussian order-from-disorder potential predicts free 
energy minima at a = 2i:n/& (in agreement with our 
simulation), at intermediate temperatures we observe in 
the finite size system a tendency of the spins to prefer 
a relative direction corresponding to a negative value of 
(cos( 6 a)) (see inset in Fig. 3(b)). This is presumably a 
result of nonlinear spin fluctuations around the classical 
ground state order, similarly to the effect of quenched dis¬ 
order [21]. We thus identify the transition temperature 
T<(Zg) as the location of the minimum of (cos( 6 a))(T) 


X (T,L)/L 2 -^> 



(| m(T,L)\)/L-^</ 2 


Z>{T)/L 



U(T)/L 


FIG. 4. (Color online) Finite size scaling of susceptibility 
x(T,L) = L 2 ~ v> Y x (£ > /L) as a function of £>/L and mag¬ 
netization (\m\(T,L)} = L~ v< ^ 2 Y m (^ < /L) as a function of 
£</L for Jth = 0.6J, Jtt = 1-0 and J = 1.22. Best data col¬ 
lapse is obtained with a BKT scaling ansatz and yields T<,>, 
o<,> and ??<.> as given in the text. 


which yields temperatures that are within error bars in 
agreement with the ones predicted from scaling. 

We note that in the critical phase that develops for 
T £ [T<, T>], the phase a behaves as a perfect, decoupled 
XY order parameter. Once the vortices bind at the BKT 
transition T>, the ensemble of thermodynamically acces¬ 
sible states divides up into distinct degenerate subspaces, 
each defined by a pair of winding numbers {n x ,n y } with 

m = / ^v 2 a(x), (l = x,y ), (7) 

where L is the linear size of the system, indicating the 
presence of an emergent topological phase [28] . The mul¬ 
tiple degeneracies of this state confirm the Polyakov hy¬ 
pothesis that a power-law phase is possible with a degen¬ 
erate vacuum. 

In conclusion, employing extensive parallel-tempering 
Monte-Carlo simulations, we have presented conclusive 
evidence for an emergent critical phase in a 2D isotropic 
classical Heisenberg spin model at finite temperatures 
thereby realizing the Polyakov conjecture [3]. Using fi¬ 
nite size scaling we have shown that the transitions are 
in the Berezinskii-Kosterlitz-Thouless universality class. 
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At low temperatures, we find direct evidence of long- 
range order in the relative orientation of the spins via 
breaking of a discrete six-fold symmetry induced by an 
order-from-disorder potential. Direct numerical analysis 
of the spin stiffness tensor, the metric of the associated 
SO(3) x U(l) topological manifold, and its Ricci flow 
will be the topic of future work. 
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Analytic calculation of the order-by-disorder potential 

In this appendix we summarize our calculation of the order-by-disorder potential terms. The aim of this calculation 
is to determine the dependence of the free-energy on the manifold of ground-state configurations (the calculation 
follows [1-3], [4]). In order to carry out the calculation we parametrize the spins as small fluctuations around one 
of the ground states. We choose the plane of spins of the triangular lattice spins in their groundstate to lie in the 
xy plane. The honeycomb spins have an order parameter, that is described by angles (3 and a. The fluctuations of 
the triangular sublattice spins are called {p, a, r) and those of the honeycomb spins are called (p,v). We choose the 
following parametrization of the fluctuations 


n A = 


1 - 


, n z 

Py+Pz 


Py 

Pz 


n B = R 


n c = R 2 


1 - 


—Z I 


Z I z 

1 _ T y +T z 


_i _ v(3 

2 2 u y 


y/3 1 _ V3 / 2 | _2\ 

n &11 4 \&y + G z ) 


"y+^a 

4 

V3(„2 , „2\ 


2 u y 4 \ J V 

O', 


( 1 I n/3- I fi+ll 

2 ' 2 'y 4 

_ l r + V3( T 2, 2) 

2 2'V ' 4 ' y ' z) 


with unit vectors 


h A = ^1 - 2 h ° + Mi/ 3 + M2« 


h B = - 


1 - 


/ip T T ^2 


\ , 

f — sin a \ 

), «= 1 

cos a 

^ 0 j 


sin /? cos a \ / cos f) cos a 

h 0 = | sin [3 sin a , f3 = cos (3 sin a 
cos (3 J \ sin (3 

and where R = R z (2tt/3) is the rotation matrix that rotates spins around the z-axis by 2tt/3. With this parametriza¬ 
tion it is a straightforward task to rewrite the Hamiltonian in the following form. The triangular spin Hamiltonian 
becomes 

H t = j tt J2 n ?■ (j2 n f + + J tt Y. n ? ■ Yy n f 

i \J(*) i'(*) / * J(i) 

here the notation j(i) denotes summation over the neighbors j of the spin with index i. This is 


H t = --NJ tt +6H t 
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with 


5H t — Jtt 

(*d> 

+J « 

<i»i> 

+ J « X! 

(*>j) 


T {piy + Piz + CT jy + a jz) ~ nPiy a jy + PizCTjz 


A (Piy Piz "h T iy “h T jz) nPiy T jy "h PizTjz 


_ f 2 , : 

4 \ a iy "+■ °»2 1 'jy 


r i7/ + Tja) 2 <J iy T iy ~h ^izF 


*2 ' J 2 


where all the linear terms have canceled out. Similarly we have 


SH h = J hh 


^ 2 L 


(m4) 2 + (m 4) 2 + (^) 2 + (^4) 2 + (^fi) 2 + (m? 2 ) 2 + (^) 2 + (4) 2 + (mS) 2 + (m? 2 ) 2 + (^) 2 + (4 


Jhh Mil 


Jhh ^ ] Mi2 


i+1+2,1 
,A _ 

z+i+2,2 


V* + ^ 




Mil ^i+i+2,1 + ^i+ 1 + 2,1 ^»+i,i 


Mi2 ^+1+2,2 + ^i+1+2,2 U i+ 1,2 


+ Mil 


+ M§ 
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,A 

i+i+2,2 ’ 


^i+2,1 ^i+1+2,1 


B C 

V i+2,2 ^i+i+2,2 



FIG. 1. Labelling of sites inside the enlarged magnetic unit-cell. 








and finally 


H th = Jth^cosapiy (pf 2 + pf_ i_ 2,2 + ^ 1 - 2 , 2 ) +cos/3sinap. jy (/4 + pf_i_ 2 A + 


— sin 
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• n / 4 r / 7 \ r \' \/3 sin <+ — cos ol (a r c \ 
+ sin pp iz (ug + vg + vg) + J th J2 - 2 - ° iy V 1 ? 2 + 4 + ^i- 2 , 2 ) 


\/3 cos a + sin a 


cos {3<Ji y (pf-[ + pg + ^ 2 , 1 ) “ sin (3ai Z {pg + pg + pf_^ t 


+ - 


— \/3 sin a + cos a 


a iy (gi+ i+ 2,2 ^*2 + u i+i ,2 


y/3 cos a + sin a 


cos/3 Viy (< 1+2,1 + 4 + 4i,l) 


,B , ,,C 


^i+i+2,1 + v il + 


l) + Jth X] “ 


a/ 3 sin a + cos a 


T *J/ {^g + i ; 2 3" 4) 


+ sin /3<jj Z ( 

i 

y/3 cos a sin ct n ( a r • n ( a r 

+- 2 - cos /3n y + ic_i,i + pgJ - sin pTi z {p a + p._ j 4 + pgJ 

-\/3 cos a — sin a 


\/3 sin a + cos a / A 

iy l ^+1+2,2 T “"i+2,2 


H- 2 - T ™ * v 

• / j, § j, , 

1 /.. , ; , A , + !/.,*, + 14 


+ ^5o + 4 ) — 


COS PTiy (uf + 1+2 l + ^ 2)1 + l/g) 


+ sin; 


i+i+ 2,1 3" ^+ 2,1 3" ^ti) • 


t-priz ( 

We transform this into fourier space with 

Pi = P{ri) = J d 2 qe lqTi p(q) 


and rewrite the resulting total Hamiltonian in matrix form. This matrix has dimensions 18 x 18. It may, however, 
be written in extremely compact block form. This will facilitate the calculation of the free energy below. The full 
Hamiltonian is 


SH = 


N 


6(2t 


f d 2 q i/> t ( JttMt JthM th \ 

J \ J th M th J hh M h J ' 


with 


M t = 


^ 1 - 3tdo 0 


0 


31- 


' J 0 J 0 


M h = 


31 

_ e -*(qi+92 )j Q 

0 

0 


_gi(9i+92)jt 

31 
0 
0 


0 
0 
31 

_ p -*( 9 1+92) 


JO 


0 

0 

_ e i(9i+92)jt 

31 


and 


M th = 


/ Jou j 0 v \ 

1 -t -t ’ 

-J6 U ~Jo v 

j 0 w 0 

V ~Jo w 0 J 
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The entries are themselves matrices: 


(l 1 1 \ 

e *(9i+®) i e iqi 

\ e *(91+®) e *g2 ^ J 


u 


V 


cos /3 


sin a 0 0 

0 sin(a + 4^) 0 

0 0 sin(a + 4/) 


— sin/3 


1 0 0 \ 
0 10 
0 0 1 / 


cos a 0 

0 cos(a + 4/) 

0 0 

The quantity if) describing the fluctuation amplitudes is given by 

1p(Ql,Q2) = (p v ,CTy,Ty,p z ,<J~,T z ,IJ.i,llf ,I/f ,ll2,ll2,V%,V2, ] '2,V2) ■ 

Finally, the correction to the free energy due to the fluctuations in a given ground state, characterized by /3 and a, is 
computed by the formula 



: ) 

cos(a + 4p) / 


5F = T 


d 2 q\og A 


with determinant 


A = det 


J tt M t J th M\ h \ 
JthMth JhhMh ) 


A straightforward analytical evaluation of this determinant is complicated by the large dimension of the matrices. 
However, we may still carry out this computation by making repeated use of the block structure of 5H. We use the 
following block matrix decomposition: 


f A B \ _ f 1 0 \ f A B\ _ f 1 0 \ f 1 B \ f A-BD^C 0 \ 

D ~ l ° 1 / v 0 I) /v 0 1 / V D ~ lc 1 ) ' 

Since all three matrices are triangular, the determinant is just 



an equality known as the Schur identity. Iterative application of this identity to the block matrices breaks the block 
structure down, until we are left with a determinant of a 3 x 3 matrix. We start with the calculation of the determinant 
for the special case /? = Here the determinant (after discarding a-independent factors) can be written as 



with J = V Jtt.Jhh and g-dependent A n . Computation of the A n shows that only n = 3 yields an a-dependence. The 
two other determinants are independent of a and are therefore not responsible for the selection of the groundstate. 
Integration over q yields the sought free energy dependence on a. Numerical evaluation of this last integral yields 


NT 

~Y~ 


d 2 q A 3 
(27r) 2 A 0 


= —10~ 4 NT cos 2 (3a) 


Jth 

~T 


6F = 


Jth 

~7 


6 
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Next we compute to lowest-order in (^j 1 ) the dependence on /?. In computing the determinant by Schur’s identity, 
one is left with a determinant of the form 


log det 


I- 


Jth 

J 


A~ 3 4 B 


Jth 


A~ Y C 


where A , B , C are 3x3 matrices. We are only interested in the lowest-order correction, which will obviously be 
proportional to (^ L ) . We find this most conveniently by rewriting the logarithm as a trace: 


log det 


Jth 


J 


-JJJ A~ l B+ -JJ a~ l c 


Jth 


J 


= Tr[log(I - 


— i A~'B 




Jth 


J 


Once the T n are known, the free energy can be expressed as 


n—1 




J 


2 n 


NT ■ 


n— 1 


Jth 


J 


2 n 


d 2 q 

W? 


T n 


Next we compute the T n by expanding the logarithm in (^j-) 


(Jf) 2n Tr (A~ 1 B+ (JfY A^CY 


Tr[log(I - 


Jth 

~7 


A~ 1 B + 


Jth 

IT 


= £(-r 


with 

Ti = Tr [A~ 1 B] 

T 2 = Tr [A^C] - ^Tr [(A^B) 2 ] 

T 3 = — Tr [A~ 1 BA~ l C] + ^Tr [(A^B) 3 ] . 

o 

Here we will only require 7\. Using the stated formula, it is straightforward to calculate the lowest-order dependence 
of the free energy on /?. We find 

SF = 0.131274 cos{2/3)NT 
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